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A new method for extracting force constants (FC) from first principles is introduced. It requires 
small supercells but very accurate forces. In principle, provided that forces are accurate enough, 
it can extract harmonic as well as anharmonic FCs up to any neighbor shell. Symmetries of the 
PCs as well as those of the lattice are used to reduce the number of parameters to be calculated. 
Results are illustrated for the case of Lennard- Jones potential where forces are exact and FCs can 
be calculated analytically, and Si in the diamond structure. The latter are compared to previously 
calculated harmonic FCs. 
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I. INTRODUCTION 

Two methodologies have been so far developed to cal- 
culate force constants (FC) of bulk crystals. One relies 
on perturbation theory and linear response formalism^, 
in which the infinitesimally small displacements of atoms 
are taken as the perturbing potential, and the required 
properties, namely the FCs, are calculated as a function 
of ground state eigenvalues and eigenstates of the crystal 
by using perturbation theory. The other more traditional 
method, called the direct method3»^i^ requires a small 
cartesian displacement of an atom placed in a supercell 
(preferably twice larger than the cutoff range of the FCs). 
Forces applied to other atoms divided by the small dis- 
placement produce the required FCs. This method, how- 
ever needs also the calculation of Born charges in the case 
of polar semiconductors, so that FCs due to long-range 
Coulombic forces can be added analytically to the ob- 
tained values. Anharmonic FCs of third or higher order 
can also be calculated in the same fashion. In the di- 
rect method, in principle two small displacements along 
a given direction are needed in order to fit the produced 
force with a second order polynomial and extract the har- 
monic and cubic coefficients. However, we are not aware 
of any calculation involving three-body (or higher) inter- 
actions. Perturbation theory, on the other hand, uses the 
so-called "2n-|-l formula" to get FCs up to order 2n+l 
from a n}^ order perturbation expansion of the electronic 
wavef unctions. One can see that the direct method, not 
being a very systematic one, soon becomes prohibitive if 
FCs of third or fourth order are needed, given their very 
large number. The calculation of fourth order, or so- 
called "quartic" terms becomes also quite involved using 
the perturbation theory method as now wavefunctions 
need to be expanded up to second order. 

In this paper, we propose a methodology to extract 



the harmonic as well as cubic and quartic anharmonic 
FCs from first-principles calculations in a systematic 
way. While harmonic FCs are used for the calculation 
of phonon spectra, vibrational modes and elastic prop- 
erties, cubic anharmonic FCs give the phonon lifetimes 
and scattering rates, and quartic ones correct the shift in 
the phonon frequencies. These quantities are the main 
ingredients in phonon transport theories which calculate 
thermal conductivity. They can also be used for con- 
structing a classical molecular dynamics potential of ab- 
initio accuracy. From such molecular dynamics simula- 
tions, thermal conductivity, which can be written as the 
heat current autocorrelation function^ and other thermo- 
dynamic properties of bulk or nanostructured materials 
can be calculated. Phase transformations can also be 
investigated by using anharmonic potentials. 

II. METHODOLOGY 

The force constants are defined as derivatives of the 
potential energy with respect to atomic displacements 
about their equilibrium position. We write the potential 
energy in the following form: 

V = Vb + ^ HjUj + ^ X! + ^ X! UiUjUk 

i ij ijk 

+ -^^^Xi]kiUiUjUkUi (1) 

ijkl 

where the roman index i labels the triplet (i?, r, a) with 
R being a translation vector of the primitive lattice, r an 
atom within the primitive unit cell, and a the cartesian 
coordinate of the atomic displacement u. In other words, 
Ui is the displacement of the atom (R, r) in the direction 
a from its equilibrium or any reference position. 4" and 
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X are respectively the harmonic, cubic and quartic FCs, 
whereas 11 is the negative of the residual force on atom i , 
and is zero if the potential energy V is expanded around 
its minimum or equilibrium configuration. In clusters or 
molecules the formalism is the same, only the translation 
vector R needs to be dropped. 

The resulting force on atom i would then be: 



= - 



dV 

dui 



ijk UjUk 



4 X! ^^J'^' UjUkUl 



jk 



(2) 



jki 



If one considers a simple FCC Bravais lattice in 3 di- 
mensions, one can discover by inspection that there are 
123 harmonic, 4867 cubic and 95663 quartic anharmonic 
FCs only within the first neighbor shell. Of course they 
are related by different symmetry properties, which will 
be used here to reduce their number to respectively to 
4,12 and 56. 



A. General symmetries of the problem 

The symmetries of the FCs are deduced from rotational 
and translational invariance of the system, in addition to 
the symmetries of the crystal itself. 

These symmetries are: 
-Invariance under the permutation of the indices: 



Xijkl — Xikjl — Xijlk — Xjikl — 



(3) 

where i, j, k and I refer to neighboring atoms. This comes 
about because the force constants are derivatives of 
the potential energy, and one can change the order of 
differentiation for such analytic function. 

-Invariance under an arbitrary translation of the sys- 
tem: 

Translational invariance of the whole system (even if not 
an ordered crystal) also implies V(u) — V(u + c) and 
F{u) = F{u + c) (which is easier to use) where u{t) are 
the dynamical variables, and c is a constant arbitrary 
shift vector (may be incommensurate with R). 

This implies the well-known "acoustic sum rule" (ASR) 
generalized to higher order FCs: 

^ V(a) => Total force on unit ceU = 

r 

= E <l.r.,H.r. yiaf3j,Run) 

= E XoXruR2r2,R,r, V {a^jS, R,R,, T^T^) (4) 



SO that diagonal terms in these tensors can be defined in 
terms of their off-diagonal elements, for arbitrary carte- 
sian components. For more details about the proof of 
these relations and the ones following on rotational in- 
variances, we refer the reader to Rcf.^. 

-Invariance under an arbitrary rotation of the system: 
Likewise if the system is rotated arbitrarily, the total 
energy and forces should not change. This leads to the 
following relations^: 

= E nor^'^e"'''': V(i^) (Torqueonunitcell = 0) 

r 

Rl,ri 



fl.2.T2 

= E Xo^tn.,n2r.,R.rAR3r3r e''^' 



\S ^■ySf j_ ff,-yP _ 



Rari 



(5) 



Here e"'^'*' is the anti-symmetric Levy-Civita symbol, 
and (Rt)" refers to the cartesian component a of the 
reference position vector of the atom t in unit cell defined 
by R. Moreover, an implicit summation over repeated 
cartesian indices is implied. 

As we see, rotational invariance relates the second to 
the third order terms, and the third to the fourth order 
terms, implying that if the expansion of the potential en- 
ergy is truncated after the fourth order terms, we need to 
start with the fourth order terms, and application of ro- 
tational invariance rules will give us constraints on third 
and second order FCs respectively. 



B. Point / space group symmetries 

The other constraints come from symmetry operations, 
such as lattice translation, rotation, mirror or any sym- 
metry operation of the space/point group of the crys- 
tal/molecule which leaves the latter invariant. 

Invariance under a translation of the system by 
any translation lattice vector R implies the following 
relations: 



ng^, = Ho", V(i?r«) 



a(3 

Rt,RiTi 



Q/3 

Qt,Ri-Rti 



^ Rt,RiTi.R2T2 ~ ^0t,Ri-Rti.R2-Rt2 
XRT,RiTiM2T2,R.:iT3 ^ X0t,Ri-R.Ti,R2~Rt2,R3~Rt3 

So in an ideal crystal, this reduces the number of force 
constants considerably (by the number of unit cells, to 
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be exact), meaning that we will use for the atoms in all 
the other cells the same FCs as those we have kept for 
the atoms in the "central" cell. 

If a rotation or mirror symmetry operation is denoted 
by 5", we must have: 

a' 
a'fS' 

^th.^Sr. = E 5^3,0, (7) 

a' p'^' 

af^jS \ ^ a' p' ^' S' C C C C 

XsT.,STuST2,Sr3 - XT,Ti,T2,T3'^a,a' Oi3,p> O^^Y Os^S' 

a' f3'^' 5' 

where Sa.a' are the 3x3 matrix elements of the symmetry 
operation 5. These symmetry relations impose a linear 
set of constraints on the force constants. 

We should note that any physically correct model of 
force constants must satisfy the invariance relations. On 
the other hand, we do approximations by truncating 
the range of FCs and their order in the Taylor expan- 
sion of the potential. Therefore imposing the constraints 
will move their value away from the true value, but has 
the advantage that they are physically correct, and will 
for instance reproduce the linear dispersion of acoustic 
phonons near fc = 0. So one should keep in mind that 
an unrealistic truncation to a too short of a range will 
produce results in disagreement with experiments. 

III. IMPLEMENTATION 

Given a crystal with its unit cell and atoms, we first 
identify its symmetry properties and construct the ma- 
trices S. Using the latter, and the equations [3] and [3 in- 
dependent FCs of each rank are identified. Then a set of 
force-displacement data calculated from first-principles 
is constructed. Next, since the data set is (better be) 
larger than the number of unknown FCs, the linear set 
of equations [H H] and [5] arc fitted with this data using a 
singular value decomposition (SVD) algorithm. In SVD 
the unknown force constants will be calculated in a least 
square sense, and furthermore linear dependencies among 
the equations will not be a problem as the result is pro- 
jected out of the "zero eigenvalue space"- . One has also 
the option to use the relations [3] and [H] to eliminate some 
of the FCs in [2] and solve for the remaining FCs. Here 
one needs to make a judicious choice of FCs to be elimi- 
nated. We prefer not to eliminate the FCs and keep the 
option of checking the violation of the translational and 
rotational invariances if only relations [2] are used. 

The data set can be obtained in several ways: a molec- 
ular dynamics run with small initial displacements, ran- 
domly moving all atoms by small displacements a few 
time steps and calculating the forces on all atoms, and 
finally displacing one atom at a time symmetrically about 



its equilibrium position by a small displacement along x,- 
x,y,-y,z and -z directions and calculating the forces on all 
other atoms. Experience has shown that the latter works 
better for the computation of harmonic force constants 
and two-body force constants in general. For three and 
four-body FCs one needs to displace at least two and 3 
atoms at a time, respectively. 

We must add that the data set is obtained not from 
a unit cell but from displacements performed in a super- 
cell whose size better exceed the range of FCs; otherwise 
the contribution of images from neighboring supercells 
might also be included in a considered FC. In some cases, 
this could lead to errors in the evaluation of FCs. No- 
tice that the exact outcome of this procedure includes, 
strictly speaking, the contribution of images as well. For 
instance, in the case of harmonic FCs, instead of $r,r' 
actually the sum $r.r' — X^l ^t,l+t', where L is a trans- 
lation vector of the supercell, will be extracted. There- 
fore supercells of low symmetry are preferable in order 
to avoid encountering FCs that can not be determined. 
For example in a cubic supercell the force constants be- 
tween the corners of the cube can never be calculated this 
way since the distance between adjacent corners never 
changes. It must, furthermore, be emphasized that the 
size of the data set as well as its accuracy are crucial in 
determining the correct set of force constants. 



IV. RESULTS 

To check the feasibility and accuracy of the method, 
we first used the Lennard- Jones (LJ) pair force, for which 
derivatives can be analytically calculated and compared 
to our results. Furthermore, the LJ forces are accurate 
within the printed number of digits by the computer, 
and do not suffer from any convergence or roundoff error 
problems. First we considered an FCC-Bravais crystal of 
LJ particles with integer cartesian coordinates. Particles 
were confined to interact only with their first neighbor. It 
was found that choosing a very small set of displacements 
(of the order of 0.001) produced the best results. Since 
in both the MD data and the fitting procedure the cutoff 
was set so that only first neighbors interact, the SVD 
procedure reproduced the exact FCs remarkably well: the 
error in the harmonic, cubic and quartic FCs was in the 
7**^, 5*^ and 4*'^ digits respectively. 

As a more stringent test, and in order to get confi- 
dence on the accuracy of the fitting method for a real 
material such as Si, we also considered a diamond struc- 
ture of LJ particles with interactions up to 10th nearest 
neighbor, but limited them in our fitting up to the 8th 
neighbor only. The energy unit was taken to be one and 
the length unit a was taken to be equal to the first neigh- 
bor distance (^3/4). No MD was performed, instead we 
only displaced the first atom by 0.0005 and 0.001 along 
the X direction and recorded the forces on all other atoms 
in a 3 X 3 X 3 cubic supercell of 216 atoms. The results are 
summarized in table ID Three-body forces are found to 
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be respectively 5 and 3 orders of magnitude smaller than 
the largest two-body values in cubic and quartic FCs. It 
is worth noticing that the fitted values obtained impos- 
ing the invariance relations are slightly worse that when 
invariances are not imposed. This is because some of the 
longer range FCs (9*'' and 10*'' neighbors) were neglected 
in the fitting procedure, and as a result the values of the 
included FCs are slightly affected when invariance rela- 
tions are imposed. In a real material case however, it 
is more important to have physically correct FCs, even 
though their values might not be exact. In any case, 
this can be checked by comparing the results with and 
without imposing invariance constraints. Imposing the 
constraints conditions in the SVD procedure resutls in 
the violation of the latter by typically 10~^, whereas the 
force-displacement data could be violated giving relative 
errors as large as a few percent, if the chosen cutoff is 
too short or first-principles forces are not very accurate, 
or displacements are too large. 

The same procedure was followed to extract FCs for 
a real material case: Silicon (diamond). In this case, 
the range of FCs is unknown and probably longer-ranged 
than we can exactly handle. The range of harmonic FCs 
was limited to 8 neighbors, that of the cubic coefficients 
was limited to 2, and quartic FCs were limited to nearest- 
neighbor interactions. First we used a data set, similar 
to the previous LJ calculation, where atoms were dis- 
placed by 0.008 and 0.016 A along all 3 directions in 
order to minimize systematic and round-off errors which 
occur in first-principles calculations. LDA-based ultra- 
soft pseudopotentials were used within the VASP den- 
sity functional simulation package^. A cutoff energy of 
300 eV and a single K-point (1/4,1/4,1/4) was used in 
the 3x3x3 cubic supercell of 216 atoms. 

This was the largest system we could handle. Per- 
haps more accurate determination of forces using more 
K-points and a higher cutoff energy would give better 
results. 

In table |TT] harmonic FCs are compared to a recent 
calculatiorti^ based on density-functional perturbation 
theoryii^. 

The phonon spectrum obtained by including harmonic 
force constants up to 8th neighbors is plotted in Fig. [1] 
versus experimental data points. 



V. DISCUSSIONS 

The test example of Lennard- Jones potential 
shows that provided the actual calculation of force- 
displacements is extremely accurate, it is possible to 
extract harmonic, cubic and quartic force constants with 
relatively good precision, by using the above method. 
Choosing low-symmetry supercells may be advantageous 
and avoids heavy calculations in big supercells. Our 
experience on Si has shown that it is possible to extract 
FCs of Silicon up to the fifth neighbor by using a small 
12 atom (1x2x3) supercell. But longer-range FCs 
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TABLE I: Comparison between analytical and numerically ex- 
tracted (SVD) values of the L J force constants in the diamond 
structure, when invariance relations|4]and[5]were not imposed, 
and only 8 neighbor shells were included in the fitting. Real 
interactions were included up to 10 neighbors {Rcut = 1.6). 
Subscripts in the first column refer to first or the second atom 
in the primitive cell (r index for $ and Rt index for 'I' and 
X ), while superscripts refer to the cartesian coordinates (a 
index) . 



would require larger cells because the fitting would fail 
in this case. In fact, large supercells are really needed for 
the extraction of longer-range harmonic FCs. Once the 
latter are known, it is possible to freeze them and extract 
higher-order FCs from a smaller size supercell but more 
accurate calculation. In this case, the force-displacement 
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TABLE II: Comparison between harmonic force constants of 
Si. Again invariance relations |4] and [5] were not imposed in 
the fitting, but were included by correcting $if . Up to 8 
neighbor shells were included in the fitting. For brevity, we 
report the results on the first 5 shells. Units are in eV/A^ 



Band Structure and DOS for Si (up to 8th neighbor FC) 
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1: Phonon band structure of Si using up to 8"* neighbor 



data can come from a molecular dynamics run of small 
amplitude, as more than one particle is needed to move 
in order to extract three and four-body terms. 

In case the so-developed potential is going to be used 
for performing large-scale molecular dynamics (MD) sim- 
ulations, one does not even need to include harmonic 



(and anharmonic) FCs beyond a few, say 4 or 5, neigh- 
bor shells. In fact most classical force fields have 
shorter range. Although MD runs might be more time- 
consuming with this method, it has the advantage of 
higher accuracy. Lattice distorsions and defects can also 
be treated using this potential, as the new force constants 
of the new atomic configuration can be obtained from the 
already developed Taylor expansion: 



^new = * 
Xnew — X 



(8) 



in which the static displacement uq is due to external 
forces, and would be obtained by solving F{uo) — 
where the force F is calculated from Eq. [51 

An additional complication arises when one is dealing 
with polar or ionic materials. In such cases Born effec- 
tive charges need to be calculated and their long-range 
effect on the FCs calculated separately using the Ewald 
summation or multipole expansion methods. Care must 
be taken since its short-range part is already included in 
the extracted FCs. So it needs to be properly subtracted 
from the Ewald sum contribution. Needless to say, this 
problem is present in all polar materials, and can not be 
bypassed. 

Using the extracted cubic force constants, one can cal- 
culate phonon lifetimes and frequency shifts. Details of 
such calculations will be published elsewhere. 



VI. CONCLUSIONS 

To summarize, we have developed a method to extract 
harmonic, cubic and quartic force constants of any crys- 
tal from first-principles force-displacement data. The 
methodology uses symmetries of the crystal to reduce 
the number of independent FCs, and can include up to 
any number of neighbor shells in principle. It requires, 
however very accurate first-principles data, in order to 
produce reliable FCs. This method paves the way for 
the development of a new generation of interatomic po- 
tentials of ab-initio accuracy. 
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